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Abstract 

We study the Stokes-Einstein (SE) and the Stokes-Einstein-Debye (SED) relations, Dt = 
kBT/6nr]R and Dr = kBT/Si^rjE? where A and are translational and rotational diffusiv- 
ity respectively, T is the temperature, ry the viscosity, Ub the Boltzmann constant and R is the 
"molecular" radius, using molecular dynamics simulations of the extended simple point charge 
model of water. We find that both the SE and SED relations break down at low temperature. To 
explore the relationship between these breakdowns and dynamical heterogeneities (DH), we also 
calculate the SE and SED relations for subsets of the 7% "fastest" and 7% "slowest" molecules. 
We find that the SE and SED relations break down in both subsets, and that the breakdowns 
occur on all scales of mobility. Thus these breakdowns appear to be generalized phenomena, in 
contrast with the view where only the most mobile molecules are the origin of the breakdown of 
the SE and SED relations, embedded in an inactive background where these relations hold. At low 
temperature, the SE and SED relations in both subsets of molecules are replaced with "fractional" 
SE and SED relations, A ~ (t/T)-^' and A ~ (t/T)-^- where ~ 0.84 < 1 and ~ 0.75 < 1. 
We also find that there is a decoupling between rotational and translational motion, and that this 
decoupling occurs in both fastest and slowest subsets of molecules. We also find that when the 
decoupling increases, upon cooling, the probability of a molecule being classified as both transla- 
tionally and rotationally fastest also increases. To study the effect of time scale for SE and SED 
breakdown and decoupling, we introduce a time-dependent version of the SE and SED relations, 
and a time-dependent function that measures the extent of decoupling. Our results suggest that 
both the decoupling and SE and SED breakdowns are originated at the time scale corresponding 
to the end of the cage regime, when diffusion starts. This is also the time scale when the DH 
are more relevant. Our work also demonstrates that selecting DH on the basis of translational or 
rotational motion more strongly biases the calculation of diffusion constants than other dynamical 
properties such as relaxation times. 

PACS numbers: 61.20.Ja, 61.20.Gy 
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I. INTRODUCTION 



At temperatures where liquids have a diffusion constant similar to that of ambient tem- 
perature water, the translational and rotational diffusion, Dt and D^. respectively, are well 
described by the Stokes-Einstein (SE) relation {l| 



Dt 



and Stokes-Einstein-Debye (SED) relation 



Dr 



QTTrjR 



knT 



(1) 



(2) 
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where T is the temperature, t] the viscosity, ks the Boltzmann constant and R is the 
"molecular" radius. These equations are derived by a combination of classical hydro- 
dynamics (Stokes' Law) and simple kinetic theory (e.g, the Einstein relation) Re- 



cently the limits ^f jhe^ SE and SED relations have been an active field of experi- 
mental 
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291 ] research. The general consensus is 



that the SE and SED relations hold for low-molecular-weight liquids for T > 1.5Tg, where 
Tg is the glass transition temperature. For T < 1.5Tg, deviations from either one or both 
of the SE and SED relations are observed. Experimentally it is found that the SE relation 
holds for many liquids in their stable and weakly supercooled regimes, but when the liquid 
is deeply supercooled it overestimates Dt relative to r] by as much as two or three orders of 
magnitude, a phenomenon usually referred to as the "breakdown" of the SE relation. The 
situation for the SED relation is more complex. Some experimental studies found agreement 
with the predicted values of the SED relation even for deeply supercooled liquids [7, 
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relation [5|, [32, 
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_3reakdown of the SED relation to the same extent as for the SE 
351]. The failure of these relations provides a clear indication of a 
fundamental change in the dynamics and relaxation of the system. Indeed, the changing 
dynamics of the liquid as it approaches the glass transition is well documented, but not yet 



fully understood 



36, 



37 



38 



39]. 



There is a large body of evidence 
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44| that, upon cooling, the liquid does 



not become a glass in a spatially homogeneous fashion 
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481. Instead the 



system is characterized by the appearance of regions in which the structural relaxation time 
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can differ by orders of magnitude from the average over the entire system 49|. The hquid 
is characterized by the presence of "dynamical heterogeneities" (DH), in which the motion 
of atoms or molecules is highly spatially correlated. The presence of these DH is argued to 
give rise to the breakdown of the SE relation 1^. Since the derivation of the Einstein 
relation assumes uncorrelated motion of particles, it is reasonable that the emergence of 
correlations could result in a failure of the SE relation. The aim of the present work is to 
assess the validity of the SE and SED relations in the SPC/E model of water, and consider 
to what extent the DH contribute to the SE and SED breakdown. 
Computer simu 
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51, 



53|, 
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ations have been particularly useful for studying DH (e.g., see Refs. [50|, 
55l |) since simulations have direct access to the details of the molecular 
motion. For water, the existence of regions of enhanced or reduced mobility has also been 
identified 
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561]. In particular, Ref. |5^ identifies the clusters of molecules with greater 



translational (or center of mass) mobility with the hypothesized "cooperatively rearranging 



regions" of the Adam-Gibbs approach 
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Most computer simulation studies on DH describe these heterogeneities based on the par- 
ticle or molecule translational degrees of freedom. We will refer to these DH as translational 
heterogeneities (TH). For water, it is also necessary to consider the rotational degrees of free- 
dom of the molecule. Recently, some computer simulation studies on molecular systems de- 

nnnnnnr 

scribed the DH based on the molecular rotational degrees of freedom [4, ll8|, l20|, l59|, l60|, iMl, l62l . 
We will refer to these DH as rotational heterogeneities (RH). For the case of a molecular 
model of water, RH were studied 62] and it was found that RH and TH are spatially corre- 
lated. This work extends those results. We find support for the idea that TH are connected 
to the failure of the SE relation, and further that RH have a similar effect on SED relation. 
Additionally, we find that the breakdown of these relations is accompanied by the decoupling 
of the translational and rotational motion. 

This work is organized as follows. In the next section we describe the water model and 
simulation details. In Sec. [Tllland Sec. [IV]we test the validity of the SE and SED relations and 
their connection with the presence of DH, respectively. The decoupling between rotational 
and translation motion is studied in Sec. |Vl In Sec. |Vl] we explore the role of time scale in 
the breakdown of the SE and SED relations and decoupling of rotational and translational 
motion. We summarize our results in Sec. IVIIi We have placed some technical aspects of 
the work in appendices to facilitate the flow of our results. 
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II. MODEL AND SIMULATION METHOD 



We perform molecular dynamics (MD) simulations of the SPC/E model of water |63 |. 
This model assumes a rigid geometry for the water molecule, with three interaction sites 
corresponding to the centers of the hydrogen (H) and oxygen (O) atoms. Each hydro- 
gen has a charge qn = 0.4238 e, and the oxygen charge is qo = —2.0 qn, where e is 
the magnitude of the electron charge. The OH distance is 1.0 A and the HOH angle is 
109.47°, corresponding to the tetrahedral angle. In addition to the Coulombic interactions, 
a Lennard- Jones interaction is present between oxygen atoms of two different molecules; the 
Lennard- Jones parameters are a = 3.166 A and e = 0.6502 kJ/mol. We use a cutoff distance 



of 2.5cr = 7.915 A for the pair interactions and the reaction field technique 64 1 is used to 
treat the long range Coulombic interactions. 

We perform simulations in the constant particle number, N, volume, V, and temperature 
NVT ensemble with N = 1728 water molecules and fixed density p = 1.0 g/cm^. The 
values of the simulated temperature are T = 210, 220,^30, 240, 250, 260, 270, 280, 290, 300, 



310, 330 and 350 K. We use the Berendsen method [65| to keep the temperature constant. 



We use periodic boundary conditions and a simulation time step of 1 fs. To ensure that 



simulations attain a steady-state equili 
least the duration specified by Ref. 



Drium, we perform equilibration simulations for at 
661]. After these equilibration runs we continue with 
production runs of equal duration during which we store the coordinates of all atoms for 
data analysis. To improve the statistics of our results, we have performed 5 independent 



simulations for each T. Ref. 



661] provides further details of the simulation protocol. 



III. BREAKDOWN OF THE SE AND SED RELATIONS 

To assess the validity of the SE and SED relations we consider a simple rearrangement 
of Eqs. ([I]) and ([2]), i.e. we define the SE ratio 

Rs. ^ ^ (3) 



and the SED ratio 



^SED = (4) 



Both RsE and -Rsed will be temperature-independent if the SE and SED relations are valid. 
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To evaluate Rse and -Rsed, we must first calculate the appropriate diffusion constants. 
Following normal procedure, we define 

where (r^(At)) is the translational mean square displacement (MSD) of the oxygen atoms 

1 ^ 

{r\At))^-J2\^,{t')-n{t)\'- (6) 

i=0 

Here, ri{t) and rj(t') are the positions of the oxygen atom of molecule i at time t and t' 
respectively, and At = t' — t. Analogously, we define the rotational diffusion coefficient 

where (y9^(At)) is the rotational mean square displacement (RMSD) for the vector rota- 
tional displacement 0i{At). Special care must be taken to calculate {ip'^{At)) so that it is 
unbounded. A detailed discussion of this procedure is provided in Appendix \M 

We also need the viscosity r] to evaluate Rse and -Rsed- Unfortunately, accurate calcula- 
tion of rj is computationally challenging. A frequently employed approximation exploits the 
fact that 7] is proportional to the shear stress relaxation time, Tg, via the infinite frequency 
shear modulus. Goo, which is nearly T-independent [67]. Additionally, we expect that 
(a "collective property") should be nearly proportional to other collective relaxation times, 
such as the relaxation time r defined from the coherent intermediate scattering function, 
F{q,At), where q is the wave vector. Therefore, we substitute rj by r, which should only 
affect the value and units of the constants in the Rse and -Rsed- For the purposes of our 
calculations, we define r by fitting F{q, At) at long times with a "stretched" exponential 

F{q,At)^exp[-iAt/r)% (8) 

where < P < 1, and we focus on the q value corresponding to the first peak in the static 
structure factor S{q). 

Now that we have the necessary quantities, we show -Rse and -Rsed in Fig. [TJ^a) and 
Fig. [11(b) with the curves labelled with "all" . Both quantities deviate at low T from the 
corresponding constant values reached at high temperature indicating a breakdown of both 
the SE and SED relations. 
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Whether there is a breakdown of the SED relation in experiments is not clear. While 
some experiments measuring dipole relaxation times show that the SED relation holds down 



to the glass transition 30|, [Sll, other experiments 



68 



show that the SED relation fails for 



low T. Our simulations are in agreement with the breakdown of the SED ratio observed 



in, e.g., Ref. [22]. Figures [T]( a) and [T](b) also show R^e and -Rsed for different subsets of 



molecules to examine the role played by DH. This is discussed in the following Section. 

IV. ROLE OF DYNAMICAL HETEROGENEITIES 
A. Identifying mobility subsets 



Many theoretical approaches (e.g. [10|, uM) attempt to explain the breakdown of SE 
and/or SED in terms of DH. To this end, we must first describe the procedure used to select 
molecules whose motion (or lack thereof) is spatially correlated. A variety of approaches 
have been used to probe the phenomenon of DH. Here we use one of the most common tech- 
niques: partitioning a system into mobility groups based on their rotational or translational 
maximum displacement. 

For the TH, we define the translational mobility, /Xj, of a molecule z at a given time to and 
for an observation time At, as the maximum displacement over the time interval [to, to + At] 
of its oxygen atom 

fii{to, At) = max{|r;(t + to) - fi{to)\ , to < t < to + At} . (9) 



For the RH, following [62l], we define a rotational mobility that is analogous to the transla- 
tional case. In analogy with Eq. ([9]), we define the rotational mobility at time to with an 
observation time At as 

iji{to, At) = max{\^i{t + to) - ^i{to)\ , to < t < to + At} . (10) 

We identify the subsets of rotationally and translationally "fastest" molecules as the 7% 
of the molecules with largest ipi and /ij, respectively. Analogously, we identify the subsets 
of rotationally and translationally "slowest" molecules as the 7% of the molecules with 
smallest ipi and /i,, respectively. The choice of 7% is made to have a direct comparison with 

6^ . but the qualitative details of our work are unaffected by modest 



the analysis of Ref. 



55 



changes in this percentage. In the following, we will refer to these subsets of molecules as 



TH and RH, fastest and slowest, depending on whether we consider the top or the bottom 
of the distribution of mobihties. We will see that comparing the fastest and the slowest 
molecules will reveal new features of DH. 

B. SE and SED relations for fastest and slowest molecules 

Having identified subsets of highly mobile or immobile molecules, we can calculate the 
ratios Rse and -Rsed by limiting the evaluation of Dt, and r to these subsets. This 
is relatively straightforward for the diffusion constants, since they depend only on single 
molecule averages. For r, the situation is more complex since F{q, At) includes cross- 
correlations between molecules. Hence we specialize the definition of F{q, At) for the TH 
and RH subsets by introducing a definition that captures the cross-correlation within subsets 
and between a subset and rest of the system. We call this function -Fsubset(Q', At), which we 
discuss in detail in Appendix [Bl 

We show the value of -Rse and -Rsed in Fig. [Dj^a) and [U^b) for the cases when only the 
fastest and slowest subsets of molecules are considered. Like the total system average, both 
the SE and SED ratios for the subsets deviate at low T from the corresponding constant 
value reached at high temperature. Therefore, we observe that the breakdowns of both the 
SE and SED relations occur not only in the subset of the fastest molecules, but also in the 
slowest. We have also confirmed a breakdown in intermediate subsets. 

The most mobile subset of molecules has a consistently greater value of -D^r/T and D^t/T 
than the rest of the system, while the ratios for the least mobile subsets are always smaller. 
This is a result of the fact that the means by which we select the different subsets most 
strongly affects the diffusion constant (see Appendix [B]), and hence the differences in the 
SE and SED ratios between the full system and the subsets are dominated by the diffusion 
constant, rather than by the relaxation time. 

In order to compare the relative deviations of these curves from the SE and SED pre- 
dictions, we normalize -Rse and -Rsed by their respective high temperature values [Fig. [T](c) 
andlU^d)]. We observe that there is a collapse of all the curves; thus, we conclude that 
both the most and least mobile molecules contribute in the same fashion to the breakdown 
of SE and SED. Moreover, this result supports the scenario that the deviation from the 
SE and SED relations cannot be attributed to only one particular subset of fastest / slowest 
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molecules, but to all scales of translational and rotational mobility. We have confirmed this 
by looking at subsets of intermediate mobility (not shown). Therefore, we interpret our 
results as a sign of a generalized breakdown in the system under study, in contrast to a pic- 
ture where only the most mobile molecules are the origin of the breakdown of SE and SED, 
embedded in an inactive background where the SE and SED equations hold (see e.g. 



These results are consistent with the results of Ref. 
via a different analysis. 

C. Fractional SE and SED relations 



22|, who arrived at the same conclusion 



When the SE and SED relations fail, it is frequently observed that they can be replaced 

3, y, y , y, y , 



by fractional functional forms 



72 



,173|,17^ 



D.-(^r\ D.^r^Y^^ (11) 



rj-i J ' ' \^ rj-i 

with < 1 and < 1. Hence we test to what degree Eqs. ffTTl) hold for our system. In 
Fig. [2] we show a parametric plot of diffusivity versus r/T for the entire system, and for the 
fastest and slowest molecules composing the TH and RH. The results at low temperature 
are well fit with the fractional form of SE and SED relations. From Fig. [21 C,t for TH is 
0.83, 0.84, 0.84 for fastest, slowest, and all, respectively, so all TH have approximately the 
same exponent. Similarly, for RH we find that is 0.75,0.76,0.75 for fastest, slowest, and 



all, respective^ 
Reference 



22 



22l | found a stronger form of this fractional relation. Specifically, Ref 
examined an "ensemble" of systems of the ST2 water model at the same T, which by 
statistical variation have fluctuations in the SE and SED ratios. Nonetheless, all systems 
collapsed to the same master curve when plotted in the parametric form shown in Fig. [2l 
meaning that the systems dominated by mobile or immobile molecules collapse to the same 
curve. While Ref. 22| employed a very different method (small systems followed for shorter 
times), the conclusion of our Fig. [2] is the same: a generalized deviation from SE and 
SED. However, Fig. [2] clearly shows that we do not find a general collapse in our present 
calculation. To understand why, we return to the fact that the method by which we define 
mobility affects much more strongly the diffusion constants than the coherent relaxation 
time, r. As a result, it is impossible to have the data for the mobile and immobile subsets 
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to collapse to a single master curve. To observe the same collapse, presumably one needs a 
more "neutral" method for selecting the mobile particles — that is one that does not explicitly 
bias toward a specific property. Unfortunately, such an approach is not obvious. However, 



we reproduced the ensemble approach of Ref. 



22j, by splitting each of our 5 simulations 



into 3 trajectories. We obtain reasonable fiuctuations that allow us to test and confirm (not 



shown) the observation of collapse of Ref. 22(]. Hence, the phenomenon of homogeneous 



breakdown of SE and SED appears to be robust for the different water models. 

V. DECOUPLING OF TRANSLATIONAL AND ROTATIONAL MOTIONS 

The SE and SED relations also imply a coupling between rotational and translational 
motion. Specifically, Eqs. ([1]) and ([2]) imply that the ratio 

^ = ^ (12) 

should remain constant as a function of temperature. Since we have already seen that the 



SE and SED ratios are not obeyed, it is likely that the ratio Dr/Dt is also violated 75|. 
However, it is also possible that D^/Dt remains constant if both D^. and Dt deviate from 
their expected behavior in the same way. 

Figure [3]^a) shows Dr/Dt as a function of temperature. As T decreases, we observe that 
Dr/D t increases, which implies that the breakdown of the SED relation is more pronounced 
than that of the SE relation. 

Experiments generally do not examine the behavior of Dr/Dt since Dr is not accessible. 



Instead, Dr is usually replaced by (r^) ^ with i = 2 19|. Here, is the relaxation time of 
the rotational correlation function 

Ce{At) = {P,{cos[p{t) ■ Pit + At)])), (13) 

where Pe{x) is the Legendre polynomial of order i, and p{t) is defined in Appendix El Fig- 
ure [3]^b) shows {Ti)^^/Dt for i = 1,2. We observe that {Ti)^^/Dt also shows a decoupling 
between rotational and translational motion. However, while Dr/Dt increases upon cool- 
ing, (tp)~^ / Dt decreases upon cooling. MD simulations using an ori/io-terphenyl (OTP) 



model 



76( 1 and the ST2 water model |2j] also find an opposite decoupling of the SE and 



SED relations depending on whether Dr or T2 is used. In the simulations of OTP, it was 
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shown that the inverse relation between Dr and T2 fails due to the caging of the rotational 
motion; this caging results in intermittent large rotations that are not accounted for by the 
Debye approximation. 

Similar to the analysis of the breakdown of the SE and SED ratios, we can test whether 
DH play a strong role in the decoupling by examining the ratio D^/Dt for the different 
mobility subsets. This is slightly complicated by the fact that we can choose mixed mobility 
subsets when calculating the ratio. Figure H] shows that the ratio Dr/Dt for all choices of 
mobility subsets approximately coincide when scaled by the high temperature behavior of 
Dr/Df. This indicates that (like the breakdown of the SE and SED relations) the decoupling 
is uniform across the subsets of mobility. 



VI. TIME SCALES FOR BREAKDOWN AND DECOUPLING 



A. Time dependent SE and SED relations 

The SE and SED relations depend on D and rj, which are defined only in the asymptotic 
limit of infinite time. In contrast, the time scale on which DH exist is finite, and generally 
shorter that the time scale on which the system becomes diffusive. As a result, making 
the connection between DH and the breakdown of SE and SED expressions is difficult. To 
address this complication, we incorporate a time dependence in the SE and SED relations, 
so that we can evaluate these relations at the time scale of the DH. This point has been 
neglected so far in the literature. To define time-dependent versions of the SE and SED 
ratios, we first define time-dependent diffusivities 

and we also define time-dependent relaxation times 

/t+At 
F{q,t')dt'. (15) 

Note that Dt{At) Dt and Dr{At) —>■ Dr in the limit At 00. The definition of 
r(At) requires some explanation: T{At) is the time integral of the intermediate scattering 
function, and T(At) will be proportional to the standard relaxation time r [Eq. ([8])] in 
the limit At — > cx3. There is a constant of proportionality resulting from the stretched 
exponential form [77|. When, instead, a DH is considered, Fsuhsetil, At) [see Eq. (IB2p ] is 
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used in the computation of r(At). We choose these definitions since, in the hmit At oo, 
they converge or are proportional to the corresponding time-independent definitions. We 
will use these time-dependent quantities to examine time-dependent generalizations of Rse 
[Eq. ©] and i?sED [Eq. O]. 

B. Breakdown time scale 

Analyzing the time-dependent ratio D(At)r(At)/T (for either rotational or translational 
motion) allows one to verify quantitatively the role of the time scale in the SE/SED ratios. 
To contrast the behavior of D(At)r(At)/T with the average over the entire system, we 
define the time dependent "breakdown" ratios as follows: 

, . _ (J(Ai)T(At)/r)p„ 

''™<^') = (D(At)r(At)/TU '''^ 
where DH refers to TH or RH. If the DH are related to the breakdown of the SE and 
SED relations, then one would expect that: (i) the 6th and 6rh ratios will show the largest 
deviations from the system average behavior at the time scale when DH are most pronounced, 
i.e. approximately at a time which we denote as t*, at which the non-Gaussian parameter 
is a maximum (see Appendix [C]) . (ii) The lower the T, the larger the peak of 6dh is (in 
agreement with the fact that the DH are more pronounced as T decreases). Figure [5](a) for 
TH and Fig. [5](b) for RH, show the behavior of 6DH(At) for the fastest subset of molecules, 
for different temperatures. Both expectations (i) and (ii) agree with Fig. [51 

From Fig. [5] we can extract the time tb when 6DH(At) is a maximum. Figure [6](a) shows 
tb for each of the four subsets: TH fastest /slowest and RH fastest /slowest. If DH play a 
significant role in the breakdown of the SE and SED relations, we would expect that the 
maximum contribution to the deviation from the SE and SED relations, occurring at tb, 
coincides roughly with the "classical" measure of the characteristic time of DH, t*. Compar- 
ison of Fig. [6](a) and Fig. Mjo) for T < 280 K shows that t* is slightly larger than for the 
slowest DH, while is shorter than tb for fastest DH. Nonetheless, tb and t* are approximately 
the same, and so the largest contribution to the SE/SED ratio is on the time scale when 
DH are most pronounced. This provides direct evidence for the idea that the appearance of 
DH is accompanied by the failure of the SE and SED ratios. 
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C. Decoupling time scales 



We next directly probe the relation between DH and the decoupling of Dr and Dt. As 
discussed above, the time scale at which the DH are observable is much smaller than the 
time scale at which the system is considered diffusive. Therefore, in analogy to the previous 
section, we incorporate a time scale in the Dr/ Dt ratio so that we can compare the decoupling 
between rotation and translation at the time scale of the DH. To this end we introduce the 

where DH refers to TH or RH. 

Figure [Tl^a) shows the results for dTH{At) for the fastest subsets of molecules. For short 
times, dTuiAt) does not depend on time and temperature, since in this initial temporal 
regime the dynamics at all temperatures is ballistic, i.e., both ((/9^(At)) and (r^(At)) are 
approximately linear with (At)^. At intermediate times dTH{A.t) develops a distinct maxi- 
mum which increases in magnitude and shifts to larger observation times as T is reduced. 
The maximum occurs at the time scale where the fastest molecules of the TH and RH 
"break their cages" and enter the corresponding diffusive regimes, see Fig. EJ^b). Therefore, 
the results of Fig. El^a) also suggest that the decoupling between rotational and transla- 
tional motion is largest at approximately the same time scale at which the DH are most 
pronounced. We note from Fig. [7](a) that dxHiAt) < 1, indicating that the decoupling of 
rotational and translational motion observed in the fastest subsets of TH is smaller than 
that from the average over the entire system. As we focus in slower subsets of TH for the 
same T, we observe that the maximum in dTuiAt) decreases at any given T. 

Figure [71(b) shows dRniAt) for the fastest subsets of molecules. Similar to the behavior of 
dxHiAt), at short times d^niAt) does not depend on time nor temperature; molecules move 
ballistically in this regime. The maxima in dnniAt) at At ^ 0.1 ps for all temperatures 
are a consequence of the librational molecular motion, enhanced in this case because we are 
selecting the fastest subset of RH. At intermediate times, we observe a broad minimum in 
dnniAt) centered at At ^ t*; this minimum becomes deeper and shifts to later times upon 
cooling, suggesting that the decoupling in the fastest subset of RH is largest at approximately 
the same time scale at which the DH are more pronounced. The fact that dTuiAt) shows 
a maximum at approximately t*, while d^niAt) shows a minimum at t* is because fastest 
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subsets of RH tend to enhance the rotational motion with respect to the translational motion, 
while the opposite situation occurs for the fastest subsets of TH. We note from Fig. Wljo) 
that djiHi^t) > 1, indicating that the decoupling of rotational and translational motion 
observed in the fastest subsets of RH is larger than that found in the average over the entire 
system. 

In short, the behavior of ^^^/(At) and dnHi^t) indicates that the emergence of DH is 
correlated to the rotation/translation decoupling, just as it does for the breakdown of the 
SE and SED relations. 

VII. SUMMARY 

In this work, we tested in the SPC/E model for water (i) the validity of the SE and SED 
equations, (ii) the decoupling of rotational and translational motion, and (iii) the relation of 
(i) and (ii) to DH. We found that at low temperatures there is a breakdown of both the SE 
and SED relations and that these relations can be replaced by fractional functional forms. 
The SE breakdown is observed in every scale of translational mobility. Similarly, the SED 
breakdown is observed in every scale of rotational mobility. Thus our results support the 
view of a generalized breakdown, instead of a view where only the most mobile molecules 
are the origin of the breakdown of the SE and SED relations, embedded in an inactive 
background where these relations hold. 

We also found that, upon cooling, there is a decoupling of translational and rotational 
motion. This decoupling is also observed in all scales of rotational and translational mobil- 
ities. In agreement with MD simulations of an OTP model |76'], we find that an opposite 
decoupling is observed depending on whether one uses the rotational diffusivity, D,., or the 
rotational relaxation time, T2. In the first case, rotational motion is enhanced upon cooling 
with respect to the translational motion, while the opposite situation holds when choosing 
T2. This is particularly relevant for experiments, where typically only T2 is accessible. 

We also found that as the decoupling of Dr/Df increases, the number of molecules be- 
longing simultaneously to both RH and TH also increases. This is counter-intuitive since a 
stronger decoupling would suggest less overlapping of TH and RH. Therefore we conclude 
that the decoupling of D^/Dt is significant even at the single molecule level. 

We also explored the role of time scales in the breakdown of the SE and SED relations 
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and decoupling. To do this we introduced time dependent versions of the SE and SED 
expressions. Our results suggest that both the decoupling and SE and SED breakdowns 
are originated at the time scale corresponding to the end of the cage regime, when diffusion 
starts. This is also the time scale at which the DH are more relevant. 

Our work also demonstrates that selecting DH on the basis of translational or rotational 
displacement more strongly biases the calculation of diffusion constants than the other dy- 
namical properties. If appropriate care is taken, this should not be problematic, but it does 
make apparent that an alternative approach to identify DH would be valuable. This is espe- 
cially true when contrasting behavior of diffusion constants and relaxation times, as is the 
case for the SE and SED relations. 
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APPENDIX A: EVALUATION OF THE ROTATIONAL MEAN SQUARE DIS- 
PLACEMENT 



To calculate [Eq. ([7j)] we consider the behavior of the normalized polarization vector 
Pi{t) for molecule i (defined as the normalized vector from the center of mass of the water 
molecule to the midpoint of the line joining the two hydrogens). The molecular rotation will 
cause a rotation of piif). A naive definition of angular displacement as Piit) — Pi(0) would be 
insensitive to full molecular rotations, since it would result in a bounded quantity. Following 



Ref. 



591 ]. we avoid this complication by defining the vector rotational displacement in the 



time interval [t, t + At] as 

iPiiAt) = / A^,{t')dt', (Al) 



where A0i{t') is a vector with direction given by Pi{t') x pi[t' + dt') and with magnitude 
given by \A0i(t')\ = cos~^(pj(t') ■ Pi{t' + dt')) , i.e., the angle spanned by Pi in the time 
interval [t',t' + dt']. Thus, the vector 0i{At) allows us to define a trajectory in a 3D space 
representing the rotational motion of molecule i, analogous to the trajectory defined by 
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fi(At) for the translationa 
displacement (RMSD) [18 



case. We define, in analogy to MSD, a rotational mean square 



59 



62| 



1 ^ 

j=0 



(A2) 



Using this form, we define Dr as given by Eq. ([71), analogous to the definition of Dt. We have 
verified that there is no qualitative difference, in the results of the present work, when the 
polarization vector is replaced by the other two principal directions of the water molecule. 



APPENDIX B: CORRELATION FUNCTIONS FOR DYNAMICAL HETERO- 
GENEITIES 

We introduce a MSD, (r^(At)), for the fastest and slowest subsets of molecules by limiting 
the sum in the Eq. (j6]) to the molecules in the corresponding subset. The different MSDs 
at T = 210 K are shown in Fig. [8](a). We note that since the most and least mobile 7% 
of the molecules will generally vary as a function of time, the molecules used to calculate 
(r^(At)) will change with time; in other words, when a molecule ceases being part of a DH, 
it is no longer considered in the computation of the MSD and the focus is shifted to the 
new subset of molecules belonging to the DH considered. Analyzing the {r^{At)) for the 
collection of subsets from most mobile to least mobile has the advantage that the mean of 
(r^(At)) over the subsets converges to the MSD for the full system. In a similar fashion the 
RMSD, ((y9^(At)), is calculated also for the fastest and slowest rotationally mobile molecules 
[Fig.Hb)]. 

To complement the single particle dynamics determined by (r^(At)) and (y9^(At)), we 
also evaluate the coherent intermediate scattering function 

^ N N 

Fiq, At) ^ -—- J2 e-'^'-^(*+^*) J] e^'^^'^W, (Bl) 
j=i k=i 

where S{q) is the structure factor. F{q, At) reflects two-particle temporal correlations in- 
stead of single-particle correlations (as in the case of the MSD). The normalization factors 
ensure that F{q, 0) = 1. In analogy to our analysis of (r^(At)), we would like to evaluate the 
contribution to F{q, At) made by subsets of molecules. Naively, one might think this can 
be simply done by limiting the sums in Eq. (1B1|) to solely those molecules within the subset. 
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However, taking the mean over the subsets of such a quantity will not recover the complete 
At), since there will be no information on the cross-correlations between the subsets. In 
order to include these correlations and define a function that, when averaged over subsets, 
will return F{q,At) (as is the case for MSD and RMSD), we simply limit one of the two 
sums to the subset, while the other sum still extends over all molecules. Mathematically, 
we define 



subset 



iVsubset Siq) ^ ^^^^^^ 
Note that one must make the choice whether to limit the sum to the subset at time t or 
t + At; we have found that in practice there is little, if any, qualitative difference in this 
choice. Thus we measure the correlations between the subset of molecules at time t with 
all molecules at time t + At. Additionally, -Fsubsetlo', 0) is not necessarily 1; forcing this 
normalization would not satisfy the desired condition that the mean over subsets returns 
the average over all molecules. In all cases, we evaluate Fsubset(Q', At) at g = 18 nm~^, the 
value of the transferred momentum at the first maximum of the structure factor where the 
relaxation is slowest (except for the g — > limit). Figure [8](c) and [8]^d) show F{q,At) for 
all molecules, and for the fastest and the slowest TH and RH. 

At this point, it is important to compare the behavior of (r^(At)) and ((y9^(At)) with that 
of F{q, At) for the TH and RH subsets. Since we define mobility on the basis of displacement, 
the behavior of (r^(At)) and {(f'^{At)) for the subsets are much more strongly affected than 
Fsuhsetil, ^t) for the subsets. Additionally, Fsuhsetil, ^t) includes cross-correlations both 
within and between subsets that a single particle definition of mobility does not include. 
More specifically, the data in Fig. [8]atT = 210K show that there is roughly two orders of 
magnitude difference between (r^(At)) for the most and least mobile molecules (and similar 
difference for (<y9^(At))). We also find that there is roughly also two orders of magnitude 
difference between the most and least mobile molecules for Dt and Dj.. For higher T, the 
difference is less pronounced. When we examine the relaxation of F{q, At) for the most 
and least mobile subsets, we find only a difference of a factor of ~ 2 between the time 
scales for relaxation. Therefore — not surprisingly — selecting mobility based on single 
particle displacement results in a much stronger effect on diffusion than it does for collective 
relaxation phenomena. This fact is important for the comparison between this work and a 
previous work 22 1. 
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APPENDIX C: CHARACTERISTIC TIME OF DYNAMICAL HETERO- 
GENEITIES 



Since we analyze the DH both in the context of translational and rotational motions, it is 
natural to ask at what time scale the TH and RH are more pronounced and to what degree 
the TH and RH subsets overlap each other. References Ssl and 6^ show that the fastest 
subsets of TH and RH form clusters, and that these clusters are larger at approximately 
the time t* corresponding to the onset of the diffusive regime, as indicated by (r^(At)) and 
(y9^(At)) respectively. Normally t* for the translational case is defined as the maximum in 



the non-Gaussian parameter 



78| 



where {r^{At)) and (r^(At)) are the fourth and second moment of the displacement distri- 
bution, respectively (the last is also the MSD). 0:2 (At) is known to be identically zero for 
a Gaussian distribution, and thus it signals when the dynamics does not generate such a 
Gaussian distribution of displacements. In the present study, we use either translational, 
filAt), or rotational, (pi{At), displacement for TH and RH, respectively, when computing 
0:2 (At). Figure E](b) shows t* as a function of T defined for the fastest and slowest subsets of 
both the TH and RH. We also include the corresponding values of t* for the entire system. 
Figure [6](b) shows that there is no qualitative difference in shape of the curve of t* (T) for 
the different subsets considered and the entire system. 

Since the values of t* for TH and RH are similar, we expect that there is some coupling 
between TH and RH. Previously, Chen et al. 79|] found that there is coupling between 
translational and rotational motion at large transferred momentum q. The maximum cor- 
relation occurs at the cage relaxation time, t*, for large values of q. Ref. [g^] found a spatial 
correlation between RH and TH. Along similar lines, we examine the overlap between these 
subsets. Figure [9] shows the overlap between the fastest subset of molecules belonging to TH 
and RH, as a function of At and T. Specifically, we count the number of fastest molecules 
belonging simultaneously to TH and RH as a function of observation time At. Similar to 
Fig. 9 in Ref. 79| , the strength of this coupling reaches its maximum at the cage relaxation 
times, but these times are consistently shorter than those reported in 7t| ; this is likely to be 
due to the fact that we consider fastest TH and fastest RH in this calculation, while Ref. [791 
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considers all the molecules of the system. Figure [9] indicates that, at the lowest temperature 
simulated, about 45% of the molecules comprising the fastest subset of TH coincide with 
the ones in the fastest subset of RH. 
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FIG. 1: (a) SE ratio, Dtr/T, for the 7% most translationally mobile molecules ("fastest"), for 
the 7% least translationally mobile molecules ("slowest"), and for the entire system (all). There 
is a breakdown of the SE relation (constant SE ratio) at low temperatures in both the fastest 
and slowest subsets, as well as for the entire system, (b) SED ratio, Drr/T, for the 7% most 
rotationally mobile molecules, for the 7% least rotationally mobile molecules, and for the entire 
system (all). Similar to (a), there is a breakdown of the SED relation (constant SED ratio), (c) 
and (d) Normalization of the curves in (a) and (b), respectively, by the corresponding quantities 
at r = 350 K. The collapse of these curves demonstrates that the relative deviations from the SE 
and SED relations are approximately the same for the corresponding mobility subsets. 
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FIG. 2: (a) Power law fits of translational diffusivity Dt as functions of t/T, Dt ~ (r/T)"^*, for 

the eight values of temperature T = 210 . . . 280 K (but not for the remaining values T = 290 . . . 

350 K), for fastest TH, slowest TH, and all molecules. We estimate S,t ~ 0.84. The dot-dashed 

line represents the normal SE behavior (^^ = 1). Consistently with the results of Fig. [H the 

deviation of these three curves from the SE behavior is almost identical as reflected in the values 

of these fractional exponents, (b) Power law fits of rotational diffusivity, D^, as functions of t/T, 

Dr ~ (r/T)"^'', of simulations in the same temperature range of (a) for fastest RH, slowest RH, 

and all molecules. We estimate ~ 0.75. The dot-dashed line represents the normal SED behavior 

iCr = !)• Also for RH, a fractional law is found with the same exponents for the three families 

considered, and, noticeably, the deviation from the normal case (^^ = 1)) is stronger for Dr than 
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FIG. 3: (a) The ratio of rotational and translational diffusivities Dr/Dt as a function of tempera- 
ture. As temperature decreases, this ratio increases indicating a decouphng between rotation and 
translational motion. The deviation of Dr is stronger than that of Dt- The line is a guide for the 
eye. (b) Same as (a) where the rotational diffusivity, Dy, is replaced by the inverse of the rota- 
tional relaxation time, with £ = 1,2, as usually done in experiments. An opposite decoupling is 
observed in (a) and (b). The lines are guides for the eye. 
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FIG. 4: The ratio of rotational and translational diffusivities, and Dt respectively, for the 
following choices of subsets: for fastest TH divided by Dt for fastest TH (0), for slowest 
TH divided by Dt for slowest TH (a), D^ for fastest RH divided by Dt for fastest RH (O)) 
for slowest RH divided by Dt for slowest RH (□). The values were normalized by the T = 350 K 
values for every curve. The fact that for these four cases D,./Dt deviates from unity (dashed line) 
to approximately the same degree indicates that the decoupling occurs across all four cases. 
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FIG. 5: (a) Time-dependent extension, bTu{^t), of the SE relation for the fastest TH at different 
T. For the sake of clarity the curve corresponding to T = 290 K was removed, (b) Time-dependent 
extension, 6RH(At), of the SED relation for the fastest RH at different T. For the sake of clarity 
the curve corresponding to T = 290 K was removed. In both (a) and (b), the maxima occur at the 
time scales corresponding to the end of the cage regime, when DH are more pronounced. These 
maxima increase upon cooling, as the DH become more pronounced. 
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FIG. 6: Temperature dependence of (a) i;,, the time at which the time-dependent extensions of 
the SE and SED relations, ^dHi have maxima, and (b) t*, the time at which the non-Gaussian 
parameter, 02 (At), reaches a maximum, t* indicates the time scale corresponding to the end of 
the cage regime. We show the results when considering molecules belonging to TH, RH, and also 
for the entire system. 
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FIG. 7: (a) Temporal behavior of the ratio of the time-dependent rotational diffusivity and trans- 
lational diffusivity for fastest TH, normalized by the average over the entire system. We show all 
the simulated temperatures, (b) Temporal behavior of the ratio of the time-dependent rotational 
diffusivity and translational diffusivity for fastest RH, normalized by the average over the entire 
system. We show ah the simulated temperatures. The figure shows that the decoupling of rotation 
from translation is increasingly more pronounced as T decreases and is a maximum (a) or minimum 
(b) on the time scale of the DH. 
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FIG. 8: Example of time correlation functions limited to subsets of DH. (a) MSD and (b) RMSD 

at T = 210 K for the fastest and slowest TH and RH respectively, as well as for the entire system. 

Intermediate scattering function, F(q,At), at T = 210 K for (c) the fastest and slowest TH, and 

entire system and (d) the fastest and slowest RH and the entire system. 
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FIG. 9: Fraction of molecules belonging simultaneously to both fastest TH and fastest RH versus 
observation time At, at different temperatures. This overlap of fastest TH and fastest RH is 
maximum at the end of the cage regime and increases upon cooling. It is almost 45% at the lowest 
T. 
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